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Abstract. We report results for reaction and vibrational quenching of the collision 
D with para-H 2 (u, j = 0) and ortho-^R, j = 1) at cold and ultracold temperatures. 
We investigate the effect of nuclear spin symmetry for barrier dominated processes 
(0 < v < 4) and for one barrierless case (v = 5). We find resonant structures for 
energies in the range corresponding to 0.01 10 K, which depend on the nuclear spin 
of H 2 , arising from contributions of specific partial waves. We discuss the implications 
on the results in this benchmark system for ultracold chemistry. 


1. Introduction 

Research with cold and ultracold molecules has witnessed an explosive growth since the 
first predictions nna and observations M of ultracold molecules in the late 1990s. 
Already, several books have been published on cold and ultracold molecules, 

and numerous reviews have been written on more specialized topics [BlEl EH DUOS], 
illustrating the continuously broadening scope of this field, as exemplified by this special 
issue. 

In addition to new proposals to produce ultracold molecules (as in Ref. ra, this 
rapidly expanding field has motivated theoretical studies of the chemical structure and 
properties of ultracold molecules, e.g., for the interaction of two KRb molecules [14] . 
The extraordinary level of control of all degrees of freedom, from internal to translational 
degrees, reached in ultracold systems has allowed the investigation of the effect of the 
nuclear spin on atom-diatom scattering, as demonstrated in Ref. [T5] with KRb+K. 

In this article, we explore such questions using a benchmark system exhibiting both 
quenching and reaction, namely H 2 +D, which can react to produce HD+H. Because of 
its small mass, this system offers the opportunity to study the contributions of low 
partial waves even at cold temperatures, as opposed to heavier systems where ultracold 
temperatures are required. This particular system is fundamental in quantum chemistry, 
and is being actively investigated in experiments capable of reaching the energy range 
explored in this work [16] . The detection of the resonance features we predict would be a 
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valuable test of the ab initio potential energy surface calculation, and their dependence 
on the nuclear symmetry provide an even stricter comparison between experiment and 
theory, the need of which is discussed in Ref. S3- In addition, H 2 +D is also relevant 
to astrophysics, especially in the astro chemistry of cold interstellar clouds jT8j , and 
even in the early universe m- Although many studies involving the scattering of H 2 
with various atoms have been performed in the ultracold regime [20], 2H [22, [23] EJj, 
the effect of the nuclear symmetry has not been studied for this reactive system at 
low temperatures. However, nuclear symmetry has been considered in studies of para- 
and ortho-H 2 , e.g., the rovibrational relaxation and energy transfer in ultracold H 2 +H 2 
collisions [25] [26, 27], and in studies of vibrational quenching of O 2 when colliding with 
He [28], or spin-changing scattering between two ground state 0 2 molecules [29] . 

In Sec. [2] we give a brief description of the theoretical and numerical tools used, as 
well as the properties of this benchmark system. We present and analyze the results in 
Sec. [3] and conclude in Sec. [4] 

2. Theoretical and computational details 

We are interested in the effect of the nuclear spin symmetry of homonuclear molecules 
in reactive scattering processes at low temperatures. To this end, we consider the 
benchmark reaction H 2 +D which involves the simplest and most fundamental diatomic 
molecule, H 2 . The nuclear part of the wave function of H 2 must be anti-symmetric 
with respect to the permutation of the two protons, and since the nuclear spin wave 
function of two spin i = 1/2 nuclei with total nuclear spin / acquires a factor (— l) 21 ^ 1 
under exchange, the rotational states of molecular hydrogen are restricted: para- 
Hydrogen corresponds to the singlet spin state with (1 = 0, Mi = 0) which allows 
only even-j rotational states, while the nuclear triplet spin state of ortho-Hydrogen 
with (I = 1, Mj = 0, ±1) allows only odd-j rotational states. 

2.1. Numerical approach 

The expression for the state-to-state cross sections, integrated over all scattering 
directions, averaged over the initial rotational states of the reactant dimer, and summed 
over the final rotational states of the product, reads 



( 1 ) 


where the generic notation n = (a vj) stands for the arrangement label “a” and quantum 
numbers (vj) of the molecular H 2 states, and k n = ^/2/IjE^j is the initial momentum 
(h = 1, atomic units are used); p a the reduced mass of the binary system atom-diatom 
in the initial arrangement (a), the initial kinetic energy is i n = E — e n , and e n are 
the rovibrational energies of the dimers. E is the (total) collision energy, J is the total 
angular momentum, and £ is angular momentum for the relative motion. The primes 
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Figure 1. Ro-vibrational levels of para-H 2 (left), and HD (right). Relevant channels 
for H 2 initially in v = 1 are shown: open channels with energies smaller than 
the absolute scattering energy E for quenching and reactions, closed channels with 
energies larger than E affecting the results, and the neglected channels above a certain 
truncation energy E max . The inset shows a zoom of the levels near v = 1, with the 
scattering energy e relative to the entrance channel (v = l,j = 0). 


indicate the corresponding quantities in the exit channel, with n' = (a 'v'j') and angular 
momentum The T-matrix T^, e , ni = 5 n ' n 8£'£ — S^, t , ni is given in terms of the S'-matrix. 

We obtain an energy dependent rate constant by multiplying the cross sections 
by the relative velocity f re i. in the initial arrangement a (D + para/ortho-H 2 ) of the 
incoming scattering channel, namely 

n(-®) n(-^) j (2) 

where t’ re i. = y / 2£’kin/ ha- The total rate constant is obtained by summing over the 
appropriate final states n'. If the arrangement remains the same (o' = a), the scattering 
process corresponds to quenching of H 2 by D to lower states, and if a' ^ a , a reaction 
occurred leading to HD + H. The corresponding rate constants are simply 


KJ?(£) = 

E 


(3) 


n'^n,a'= 

a 


K*(S) = 

E 

K n t^ n (E) . 

(4) 


n'^n,a'^a 


We note that if n' = n, we obtain an elastic collision: this is discussed in more detail in 
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To solve for the S'-matrix and thus Eq.flT]), we modified the ABC reactive scattering 
code developed by Manolopoulos and coworkers [30], which is based on Delves [31] 
hyperspherical coordinates, with a separate set of Delves hyperspherical coordinates for 
each arrangement sharing the same hyperradius coordinate. The full wave function 
is expanded in eigenbases corresponding to internal motion coordinates for each 
arrangement (after the combined multiple-arrangement basis is orthogonalized), and 
the resulting hyperradial coupled-channel equations are solved using the log derivative 
method [32]. The propagation starts at p — p min (where the potential is highly 
repulsive), and ends at a sufficiently large p = p max , where the S-matrix is extracted 
by imposing asymptotic boundary conditions. We note that the ABC code is well 
suited for system involving barriers; it was tested at high energies for a number of 
benchmark systems such as H + H 2 , F + H 2 , Cl + H 2 , and their isotopic counterparts 
[331E2 S3 [361 E7j, 

and modified versions to study certain benchmark chemical reactions 
in the ultracold regime [211 SHI E2J ESI 140] . 

To obtain fully converged results is expensive computationally, and a sufficient 
number of terms in the sum in Eq. (JT]) are required for higher energies. In Sec. 3 our 
results will show clearly that a maximum value of J max = 4 is required for energies in the 
kelvin regime, i.e., for E < 10 K. For much higher energies, which are outside the scope 
of this work, the value of J max increases steadily. We remark that although the blocks 
for each J are separate (due to conservation of total angular momentum), the coupled 
problem for a given J becomes rapidly prohibitive as J increases. As described in |41j . 
we modified the ABC code’s structure and adapted it to the cold and ultracold regime; 
most important in our implementation is the ability to follow the progress of convergence 
in great detail with a save-and-restart feature, which allows us to go back and extend 
the radial propagation of a given run if needed. We typically monitor the convergence 
simultaneously for a large number of collision energies and for many initial states in a 
single run. Typically, we propagate to p max > 40 a.u. to obtain converged results using 
small integration step size (A p < 0.002 a.u.) |41j . Finally, in addition to truncating 
the sum over J (and partial waves), we also restrict the number of channels by fixing 
a maximum energy F max . As illustrated in Fig. [T] for the case of the (v = 1 ,j = 0) 
initial channel of H 2 , closed channels above A max are neglected. This figure also depicts 
the difference in the number of molecular states between para-H 2 and HD; The density 
of states in HD is larger not only because it is more massive, but also because the 
restriction on rotational states j is lifted. 

2.2. Interactions and Potential Energy Surface (PES) 

The benchmark system we consider here, H 2 +D, has already been studied at ultracold 
M and higher temperatures 02! • Accurate ab initio potential energy surfaces (PES) 
exist [43, 02] for this system. As in our previous studies SD m, we adopted the 
electronic ground state PES of Ref. [03]. 

This system possesses an activation barrier, which is illustrated by a contour plot 
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Figure 2. (a) Contour plot of the PES for H 3 with the vibrational levels v = 0, • • •, 4 
plotted for two geometries (9 = 180° and 90° as defined in the sketch above). This 
plot shows that the barrier is always present for v = 0, and partially for v = 1 and 
2. For higher v its influence is reduced, (b) Progression of the barrier as 9 varies: 
it appears between 110° and 115° for v = 1, and only near 90° for v = 2. The line 
passing through the middle shows the reaction coordinate £. (c) Energy surface along 
the reaction coordinate £ for different angle 9. The position of the vibrational levels 
v = 0 and 1 (with j = 0) for H 2 and HD are also shown; v = 0 of both is always below 
the top of the barrier, indicating that reactions require tunneling, while v = 1 of H 2 is 
above the barrier for 9 > 112°. 


of the vibrational levels in Fig. [2](a,b); basically, the barrier is present for any geometry 
for v = 0 , requiring tunneling for reactions to occur, while v — 1 is partially affected by 
the barrier (he. the barrier is present for a finite range of angles only). For higher v, 
the barrier does not block reactions (although it still affects the scattering rates). This 
figure also shows the reaction coordinates £. Following the line £, we plot the barrier 
for different angles 6 and the lowest energy levels of H 2 and HD in Fig. [2](c). We can 
clearly see that v = 0 requires tunneling to produce HD, while v — 1 is higher than the 
barrier at angles starting at 112 °. 

Hyperhne interactions are not considered in our treatment. Several studies have 
pointed out the effect of hyperhne coupling in cold and ultracold molecular dynamics, 
mainly in the inelastic processes of non-reactive systems such as YbF+He [45], NH+Mg 
[46], MnH+He [47], or in reactive systems such as Na 2 +Na in which only the atomic 
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hyperfine state is modified while the molecular hyperhne state remains unchanged [48]. 
Other studies involving H 2 at higher temperatures for the interpretation of astrophysical 
spectra, such as CN or HCN interacting with para-H 2 ra, or HC1 with para- or ortho- 
H 2 [5Q], give approximate methods to account for hyperhne interactions. However, like 
in our case, the hyperhne structure of H 2 is not considered. The H 2 hyperhne splittings 
in its electronic ground state X 1 S+ are negligible, so that the H 2 +D scattering will be 
dominated by the hyperhne structure of D (327 MHz or 15.7 mK). Since the focus of 
this work is the structure found in the range of 100 mK to roughly 10 K, the effect of 
hyperhne interaction can safely be neglected in our study of nuclear symmetry effects. 
However, at sub-millikeIvin temperatures, the hyperhne effects will become relevant, in 
particular if a resonance is located in the regime of extremely low energies. 

3. Results and discussion 

We give here an overview of the influence of the nuclear spin state on the scattering 
properties (reaction, quenching, and elastic processes) of H 2 with D. The nuclear singlet 
spin state of para-Hydrogen with (1 = 0, Mj = 0 ) allows only even-j rotational states, 
while the nuclear triplet spin state of ortho-Hydrogen with (I = 1, Mj = 0, ±1) allows 
only odd-j rotational states. We consider here only the lowest rotational states j = 0 
and j = 1 for the initial para- and ortho-H 2 , respectively. In what follows, we show 
results for scattering energies ranging from 1 pK to 100 K; as mentioned in the previous 
section, hyperhne interactions will become relevant for energies below roughly 15 mK, 
and thus our results below that energy are approximate. 

3.1. Overview 

We computed the energy-dependent rate constant K^(i7) for quenching and K^(i7) for 
reaction given by Eqs. ([3]) and (j4j) , respectively, for scattering of D on initially prepared 
para-H 2 (u,j = 0) and ortho-H 2 (u,j = 1). Here, n = (a = l,v,j = 0) for para-H 2 
and (a = 1 ,v,j = 1) for ortho-H 2 , with a = 1 labeling the arrangement H 2 +D (a = 2 
stands for the arrangement HD+H). Figure [3] shows the overview of the results for the 
six lowest vibrational levels v of H 2 . Note that, for v = 0, only the reaction channel is 
opened for both para- and ortho-H 2 . Overall, we find that K V (E) of para- and ortho-H 2 
share many features; for example, the position of resonances occur at roughly the same 
energy, since they are due to shape resonances in the entrance channel. However, we 
also notice differences, such as the values of K V (E) and the absence of certain resonant 
features in a few cases. 

We find structures in the range of scattering energies corresponding roughly to 100 
mK and 10 K. These features are due to low partial waves, with the more pronounced 
ones found at very low energies, near 15 mK, arising from p- wave resonances. At higher 
energies (above 10 K) no structure or resonance feature can be found, because the much 
higher centrifugal barrier (for higher partial waves, i > 3) in the entrance channel lifts 
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Figure 3. Rate constant K V (E) (v < 5) vs. collision energy corresponding to the range 
of 1 [i K to 100 K for reaction (top panels) and quenching (bottom panels) processes 
for para-H 2 (left panels) and ortho-EE (right panels) interacting with D. 


the shallow attractive well entirely (see Fig. [4] for D+para-hE in v — j — 0), and thus 
removes any possible bound states (i.e., no van der Waals complexes exist for large 
t). For collision energies between 10 K and 100 K, the rate constant follows a simple 
behavior (nearly constant). For much higher energies (above 100 K, not shown here) an 
exponential increase is typical (especially for lower vibrational initial states). Finally, 
we also notice that both reaction and quenching rate constants for ortho-EE are slightly 
larger than the corresponding rates for para-H 2 . This is to be expected, since internal 
energy of ortho-hE initially in j = 1 is slightly higher than that of para-kE in j = 0, 
thus reducing the effect of the reaction barrier, and hence increasing the reaction and 
quenching processes. 
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Figure 4. Effective potential for D + para-H^?; = 0 ,j = 0) for different partial 
waves t. For l = 0 and 1, the van der Waals complex of H 2 +D support a bound state 
(indicated as a straight horizontal lines in both cases). As t increases, the centrifugal 
interaction lifts the well as to prevent bound states and even resonances. 

In the next sections, we describe the results for both para- and ortho-H 2 in more 

details, but we first make a general remark about the magnitude of the rate coefficients 

in Fig. |3] As we discussed in our previous work HU, tunneling through the reaction 
barrier for v = 0 and v — 1 is responsible for the small values of K R (if) in the energy 
range considered here; in particular, for v = 1, we have K R < K c h However, for v = 2 
and higher, tunneling is no longer important; moreover, due to the restriction on the 
rotational levels for para- and ortho-H 2 , and also because of the increased reduced mass 
of HD, there are many more reaction channels ( v',j') for the product HD than for 
quenching. Consequently, we have K R > K f - for v > 2. 

3.2. Para-H 2 

In the case of para-H 2 with j = 0, we have i = J only, which reduces the summation 
over £ to a single term. The cross section (JT]) becomes 

00 \J+j' | 

<W£) = pr £ t 2 - 7 +1) £ 

K n J =0 

The resulting rate constant K V (E) = v ie ia for both reaction and quenching processes are 
shown in Fig. 0 For scattering energies ranging from 1 /iK to roughly 100 K, only J = 0 
to 4 contribute significantly to the cross sections and rate constants: these individual 
contributions are also shown. We now discuss each individual initial vibrational level v 
for both reaction and quenching. 

For v = 0, only reactions into HD(u' = 0, j' < 2) are possible for this range of 
energies: quenching to lower states is not possible since H 2 is already in its lowest state, 


Tn't'nj(E) 


(5) 
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Figure 5. Rate constant K V (E) (v < 5) vs. collision energy corresponding to the range 
of 1 [x K to 100 K for reaction (left panels) and quenching (right panels) processes for 
para-H 2 + D. Each panels shows the total rate constant as well as the contributions 
of individual J = 0 ... 4. 



















































Para/ortho-H 2 +D reactions 


10 


and reactions to higher v' or even higher j' would correspond to excitations and require 
larger scattering energies; e.g., about 500 K to excite H 2 (u' = 0 ,f = 2) and 350 K 
to excite HD(V = 0 ,j' = 3). Due to the presence of a barrier (see Fig. [2j) , reactions 
take place through tunneling, leading to a very small rate constant. We find a large 
resonant feature near 200 mK arising from the J = 1 contribution (or p- wave i = 1 
since i = J for para-H 2 initially in j = 0). As shown in Fig. [4l the resonance occurs 
because the centrifugal interaction pushes the bound state of the van der Waals complex 
near the scattering threshold. Below 10 mK, the main contribution is from J = 0 (or 
s-wave i = 0), while J = 2 (or d -wave i = 2) impacts K(E') at roughly 3-4 K. At 
higher energies, the other partial waves add up (including possible excitations), but no 
structure is present. 

For v > 0, both reaction and quenching processes have qualitatively the same 
behaviors, with similar contributions from the different J. Except for v = 1, the reaction 
rate constant is larger than the quenching rate constant, simply reflecting the larger 
number of exit channels to form HD as compared to the number of channels to remain 
as H 2 , roughly 2-4 in that range of levels v (see Fig. [[]). Level v — 1 is slightly different 
because the presence of the reaction barrier is still relevant for a sizable range of relative 
orientations: as shown in Fig. [2], the reaction barrier disappears for angle larger than 
112°. The presence of the barrier for angle less than 112° implies a reduced reaction 
rate, since it requires tunneling to occur. 

The J = 0 contribution, which corresponds to s-wave (l — J — 0) scattering for 
H 2 in initial states j = 0, dominates the rate constant at very low energies, as is well 
known. However, note that for energies above 100 mK, the s-wave contribution has a 
minimum; thus, higher J-contributions do become dominant as the scattering energy 
increases. It begins with J = 1 corresponding to p -wave scattering (l — J — 1), followed 
by J = 2 (d -wave scattering with i = J = 2), and so on. The exact energy at which 
these higher J contributions become important depends on the energy surface and the 
position of levels. 

For both v = 0 (reaction only) and v — 1, the main feature is a scattering resonance 
due to J = 1 (p-wave scattering with l — J — 1), arising from a quasi-bound state in the 
van der Waals H 2 • • -D complex. From v — 0 to v — 1, the quasi-bound level gets closer 
to the threshold, leading to a shift of the resonance to a lower scattering energy (from 
200 mK for v = 0 to 15 mK for v — 1). This quasi-bound level becomes bound for v = 2, 
removing the resonant feature from the rate constant, and becomes more deeply bound 
as v increases, moving the J = 1 contribution to higher scattering energies. The J = 2 
contribution becomes important at higher energies in the range of a few K, and shifts 
slightly to lower scattering energies as v grows; it starts to be important for v = 2 , 
and becomes the key feature at higher v. Similarly, J = 3 (/-wave scattering with 
i = J = 3) starts contributing at roughly 10 K, with the exact position of its maximum 
contribution shifting slightly with v; it leads to a perceptible feature only for the highest 
v, its contribution being otherwise a simple addition to all higher J giving the total rate 
constant. This is illustrated by J = 4 (g-wave scattering with i = J = 4) which simply 
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Figure 6. Rate constant K V (E) (v < 5) vs. collision energy corresponding to the range 
of 1 [x K to 100 K for reaction (left panels) and quenching (right panels) processes for 
ortho-H 2 + D. Each panels shows the total rate constant as well as the contributions 
of individual J = 0 ... 4. 


adds up to the total rate constant without leading to features. 
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3.3. Ortho-H 2 

In the case of ortho-H 2 with j — 1, since £ = | J — j\, • • ■, J + j, three partial waves 
contribute to a given J, except for J = 0. In fact, for J = 0, we have only the £ = j = 1 
p-wave scattering contribution, while J = 1 includes the £ = 0,1 and 2 (s, p, and d- 
waves) contributions, J = 2 the £ = 1,2, and 3 (p, d and /-waves) contributions, and 
so on. The cross section (jT]) becomes 



The resulting rate constant K V (E) = v re \a for both reaction and quenching processes 
are shown in Fig. El Many of the features and behaviors described for para-H 2 apply 
for ortho-H 2 as well; only J = 0 to 4 (also individually shown) contribute significantly 
to the cross sections and rate constants for scattering energies ranging from 1 /iK to 
roughly 100 K (recall that the hyperfme interactions omitted here will become relevant 
below roughly 15 mK). For the same reason as for para-H 2 , v = 0 leads to reaction 
only, since quenching to lower states of H 2 is not possible (para-ortho mixing due to 
nuclear spin-rotation coupling is extremely weak and is neglected here), and the value is 
very small because the reaction must occur via tunneling through a barrier. Similarly, 
except for v — 1 for which the reaction barrier still influences the rate, the reaction rate 
constant is larger than the quenching constant for v > 1 . 

However, the exact features are different from those of para-H 2 . This can be 
understood in terms of the different couplings between the various states. We now 
describe in more detail the results for various initial v. Except for J = 0, all other J‘ 
include contributions from three initial partial waves t. However, even if a given partial 
wave £ participates in various J-contributions, their effect is slightly different since each 
total angular moment J leads to different coupling strengths. For example, the p- wave 
£ = 1 initial partial wave occurs in the contributions of J = 0, 1, and 2. For v = 0 
(reaction only), the p- wave resonance appears only in the J = 2 contribution, while the 
p- wave scattering corresponds to a bound state in J = 0, and has almost no influence 
over the s-wave dominated J = 1 contribution. We note that the higher J are usually 
favored because of the 2J + 1 factor appearing in the expression for the cross section 
(! 6 l) . Still for v — 0, the s-wave scattering is present only in the J = 1 contribution, 
while the higher J = 3 and 4 contain £ = 2,3, and 4, and 3, 4 and 5, respectively. As 
in para-H 2 , the effect of partial waves higher than d- wave {£ = 2) is hidden in the total 
rate; they do not produce noticeable structures. 

For all other v, the J = 2 does not contain a p- wave resonance; traces of it appears 
only in the J = 1 contribution of v = 1 and 2. In more detail, for both v = 1 and 2, 
J = 0 includes only p-wave scattering without resonance. Also, J = 1 contains the s, 
p, and d- wave scattering; the position of the maximum for p- wave is shifted relative to 
the other J, illustrating the sensitivity of the position of the bound level in the van der 
Waals complex to the couplings. This is less so for the other partial waves since the 
bound states are not close to the threshold. The J = 2 contribution includes results 
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from p, d, and /-scattering, with the £ = 1 and 3 components being dominant, while the 
/ = 2 is obscured by them. This occurs for all v. As for v — 0, the higher contributions 
J = 3 and 4 are dominated by £ = 2 and £ = 3, respectively, with signals from £ = 3 
and higher being hidden in the total rate constant. For v > 2, the p-wave feature of 
J = 1 is absent, preventing a double-peak structure in the total rates. 


3-4- Elastic cross sections 


Here, we describe very briefly the results for the elastic cross section, which reads 

'2J+1\ g g 


<(E) = o*JE) = 


K V 2 J + 1 


£'=\j-j\ e=\j-j\ 


Set - Sf 


ni' n£ 


(E) 


(7) 


For para-H 2 initially in j = 0 (with £ = J) we have 


y (2j +1) 


U 2 

J=0 




and for ortho-H 2 initially in j = 1 ( with £ — \ J — 1|, • • •, J + 1) 


( 8 ) 
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J+3 

E 
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(9) 


The corresponding results are shown in Fig. [7] for v < 5. As before, there is some 
structure in the range of roughly 10 rnK to about 10 K. The cross sections for both para- 
and ortho-H 2 +D have the same order of magnitude for all v, except in cases for which 
the p -wave scattering leads to a sharp feature. This occurs for v = 0 and 1 of both 
para- and ortho-H 2 , and for v = 2 of’ ortho-H 2 . Overall, the elastic cross sections have 
basically the same value around 2000 Oq, with very small variation about that value. 
This is expected, since the entrance channels are very similar for all those initial states 
(see Fig. ISJ) , as opposed to the reaction and quenching processes where the position of 
the exit channels and the coupling strengths dictate the large range of constant rate 
values. 


3.5. State-to-state results 

In this section, we give an overview of the richness of the state-to-state processes. We 
first discuss the branching ratios of product formation from each initial vibrational 
state v of H 2 into individual vibrational levels v' (summed over j') for quenching and 
reaction separately, together with the branching ratio into all product channels (i.e., 
quenching and reaction combined). This is followed by an example of a typical case for 
rotationally-resolved rate constants. In all cases, we contrast the results of para- and 
ortho-H 2 . 

3.5.1. Products branching ratios. We first give the branching ratios for para-H 2 initially 
in the (v,j = 0) state. They are depicted in Fig. [9] the left panels show ratios for 
quenching alone (dashed lines), the middle panels for reaction alone (solid lines), and 
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Figure 7. Elastic cross sections a^ (E) (v < 5) vs. collision energy corresponding to 
the range of 1 /jK to 100 K for D + para-H 2 (left panels) and ortho-IE (right panels). 
Each panels shows the total rate constant as well as the contributions of individual 
J = 0 ... 4. 
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Figure 8. Effective potential for i = 0 for D + H 2 (v,j = 0). The asymptotes of the 
curves for the different levels v are all set to zero. The curves are very similar at large 
separation, and start to become different at short range for higher v. The difference 
becomes sizable at energies corresponding to roughly 200 K. 


the right panels for the combined inelastic processes (quenching: dashed lines, reaction: 
solid lines). We omit the initial v = 0 panels since there is no quenching possible and 
all reactions end up in v' — 0 of HD, and the panel for v — 1 of quenching alone, since 
all quenching goes into v' — 0 of H 2 in this particular case. 

We find that the branching ratios are basically insensitive to the scattering energy 
within the 1 rnK to 10 K range, even though there is structure in the corresponding rate 
constant shown in Fig. [5} This simply signifies that all the flux into the various product 
channels are changing in unison with the appearance of scattering features; these are 
due to shape resonances and higher partial waves in the entrance channel that raise all 
scattering amplitudes in the same proportion. In detail, starting from v — 5, the largest 
fraction (32%) of quenching is in the uppermost level v' = 4, followed by 28% in v' = 3, 
19% in v' = 2, 12% in v' = 1, and 9% in v' = 0. For reaction alone, v' — 5 dominates 
the reaction products formed with 41%, followed by 25% in v 1 = 4, 14% in v' = 3, 
while each of v' — 2, 1, and 0 have between 6 % and 8 %. When all product channels 
are combined, it becomes clear that the most probable products stem from the reaction 
forming HD in v' — 5 (29%), and 4 (17%), with the rest of reaction and quenching being 
comparable at 10% or less. As discussed in Section 13.61 below, this is to be expected, 
since the number of possible exit channels is larger in the reaction arrangement than in 
the quenching arrangement | |4T| . 

The same general behavior takes place for v > 1: for v — 4, both v' = 3 and 2 have 
the same quenching ratios (~32%), while v' = 4 leads the reaction ratios with 50%, and 
the most probable products are HD in v' = 4 (39%) and 3 (19%), the remaining being 
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distributed over the other channels. For v' = 3, quenching is largest into v' — 2 (39%), 
reaction into v' — 3 (60%) and 2 (22%), and the most probable products are HD in 
v' = 3 (51%) and 2 (18%). Similarly for v = 2, 62% of the quenching ends up in v' = 1 
and 38% in 0, while 64% of reaction goes into v' = 2, 23% in 1, and 13% in 0: most of 
the products are HD in v' — 2 (52%) and 1 (19%), with the rest distributed over the 
remaining channels. 

Now, the behavior for v = 1 reflects the considerable effect of the barrier on the 
scattering: 80% of the reaction into HD produces the lowest level v' = 0 as compared 
to 20% in the highest open level v' = 1. The reaction fractions are much smaller than 
the quenching fraction; specifically, 74% for quenching into H 2 (u' = 0), 20% and 6 % 
into v' — 0 and 1 of HD, respectively. Finally, for v = 0 (not shown in the figure), only 
reaction into v' = 0 of HD is possible. 

In general, we find that the products with the largest (although not dominant) 
branching ratios are into highest v' allowed: v' — v for reaction, and v' = v — 1 for 
quenching. Also, reaction products are more likely than quenching, which is expected 
since there are more possible exit channels available in the HD reaction arrangement as 
compared to the H 2 quenching arrangement. These conclusions apply to all initial states 
v, except when tunneling through the barrier is important, as discussed in Section [3761 

The same overall conclusions apply to the case of ortho-H 2 initially in (v,j = 1) 
shown Fig. (TH in general, the branching ratios are rather insensitive to the scattering 
energy in the range of 1 mK to 10 K, the largest branching ratios are into products with 
the highest v' allowed (with reaction to form HD more probable than quenching of H 2 ), 
and tunneling through the barrier affects the lowest initial levels v = 1 and 0. However, 
we note that the branching ratios are changing for both v = 2 and 1. In v — 2, the 
quenching into v' = 1 increases from 60% to 85% near 50 mK, while quenching into 
0 decreases accordingly from 40% to 15%. For the reaction products, the branching 
ratios into v' = 2 and 0 decrease slightly at the same energy, while it increases into 
v' = 1. The net effect in the combined branching ratios shows a sizable increase of the 
quenching into v' = 1 (reaching 25%) while reaction into v' = 2 decreases to 40%. For 
the initial v = 1 state of H 2 , the net result is an increase of the branching ratio into 
H 2 (w' = 0) from 60% to 75% near 200 mK, while the reaction into HD in v' = 0 and 1 are 
both decreasing accordingly. This is an example showing that para- and ortho-H 2 have 
different behaviors. To better understand the origin of those variations, we investigate 
rotationally-resolved rate constants for a specific case, namely v = 2 , in the next section. 

3.5.2. Rotationally-resolved rate constants for f / 2 in v = 2. Because of the large 
number of exit channels for a given entrance channel, we restrict this analysis to a 
typical entrance channel showing structure, while still having a manageable amount of 
exit channels. We selected to show the state-to-state inelastic (reaction and quenching) 
rate constant for para- and ortho-H 2 initially in v — 2. Other entrance channels would 
exhibit the same type of behavior, though the exact details will vary. However, as will 
become obvious below, the rapidly growing number of channels make it difficult to show 
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Figure 9. Branching ratios of product formation for D + para-EQ (v,j = 0) into 
individual vibrational levels v' (summed over j') for quenching of H 2 only (left panels: 
dashed lines) and reaction forming HD only (middle panels: solid lines). The right 
panels show the branching ratios when both quenching and reaction are combined 
(quenching: dashed lines, reaction: solid lines). The initial v = 0 panels and the panel 
for v = 1 of quenching alone are omitted (see text). 


those results. 

We begin with para-H 2 initially in (v = 2 ,j = 0), and show the state-to-state 
rate constant for quenching in Fig. [HI and for reaction in Fig. [12l For quenching, we 
show the results in v' = 1 (top left panel) and v' = 0 (bottom left panel): in each of 
exits levels, we show the individual contributions of rotational states j' to the total rate 
constant (with the quantum numbers labels following the same ordering as the partial 
rate constants). To the right of those two panels, we show the contribution of individual 
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Figure 10. Same as Fig. [9] for D + ortho-H 2 (u,j = 1). 


J total angular momenta for each channel j'. In all cases, all axes have the same range. 
For v' = 1, the most important contribution comes from j' = 2, followed by 4, 0, 6, 
and 8: each have a very similar energy dependence. The individual contribution of 
each J are shown in the right panel: they all have the same energy dependence, with 
the main contribution arising from J = 1 and 2 (at higher energies). The results for 
v' = 0 are very similar, although the ordering of the j' contributions changes a bit and 
an additional j' = 10 appears. 

For reaction (see Fig. fl2li . very similar behaviors are found, although the number 
of exit channels v' now include 2 as well, and both odd and even j' are allowed. The 
exact ordering of contributions differs from the quenching case, but the overall system 
exhibit the same dependence. 
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Figure 11. Rate constant K (E) for para-H 2 (v = 2) for quenching into specific exit 
channels v' (left panels). Each shows the contributions of individual rotational states 
j'. For each value of j', the contributions of each J = 0 ... 4 are shown (right panels). 


The same information for the case of ortho-H 2 initially in (v = 2, j = 1) is shown in 
Fig. [13] for quenching and Fig. [14] for reaction. For quenching alone, the left panels 
in Fig. [13] depict the rate constant for v' = 1 and 0, each for the various odd j' 
contributions to the total value. Although these contributions show similar dependence 
on the scattering energy for both v', one can notice that specific j' have larger increases 
than others near the resonant feature at 50 rnK. For example, in v' — 1, j' — 7 becomes 
dominant at the resonance, even though it is the smallest contribution away from it. The 
same is true to a lesser extent for the other j'. This is more apparent in v' = 0, where 
not only j' = 7 (which becomes dominant) but also 9 and 11 have sharper increases at 
the resonant structure. The origin of those variations is due to the different couplings 
between the various states, as discussed in Section 1373] Basically, this feature arises from 
a resonance in the p -wave {£ = 1) initial partial wave, which occurs in the contributions 
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Figure 12. Same as Fig. |TT] for reaction. The quantum numbers v' and j' label the 
rovibrational channels of the product HD. 


of J = 0, 1, and 2, and their effect is slightly different since each total angular moment 
J leads to different coupling strengths. 

For reaction alone, shown in Fig. [T4l we have a very similar story, though there 
are many more possible product states (including v' = 2, and even and odd j'). Again, 
the structure is due to the p-wave in the entrance channel, and the contributions of the 
various j 1 increase different near the resonant feature at 50 mK. For example, for v' = 2, 
j' = 3 and 4 have sharper increases, while the same is true for j' = 2, 6 and 8 of v' — 1. 
These “unequal” changes in the contributions of the total rate constant in each product 
v' explains the difference between the state-to-state scattering of para- and ortho-H 2 - 

Rotational product distributions for a given scattering energy E , i.e., the 
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Energy (kelvin) 


Figure 13. Rate constant K (E) for ortho-H 2 (v = 2) for quenching into specific exit 
channels v' (left panels). Each shows the contributions of individual rotational states 
j'. For each value of j', the contributions of each J = 0 ... 4 are shown (right panels). 


probabilities p™, v r(j',E) to form a product in a rotational state j' (in vibrational state 
v') in arrangement a' can be obtained by simply dividing the rate constant for this 
specific channel by the total rate constant for the corresponding product channel v r , 
namely p n a , y (j',E) = K al yj>^ n (E)/J2j'K',v',j'^ n (E), where n = (a,v,j) stands for 
the entrance channel (v,j) in the reactant arrangement a. By summing over all open 
exit channel within one arrangement a', one can find the rovibrational distribution 
for quenching (Q: a' = a) or reaction (R: a' 7 ^ a), or simply pQ^(v',j',E) = 
^Q{R),v',j'<-n(E)/KQ( r )(E). Finally, the distribution over all possible channels can be 
obtained by dividing K q{r) iV 'j’^ n (E) by total rate constant K n (E) = K ®(E) + K^(E). 
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3.6. Statistical model 

Scattering problems with rearrangement are computationally demanding when a full 
quantum mechanical approach is used; thus, approximate models have been developed, 
often based on statistical formalism and semi-classical treatments [HU 152]. Here, we 
compare the results of our full quantum computation with a statistical treatment based 
on the phase-space theory (PST) [HI] . 

Because the entrance channel plays a special role at low scattering energies, it is 
advantageous to interchange the order of the J and i summations in Eq. (JTJ), and rewrite 
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Following Miller 
(i.e., n 7 ^ n') as 


£=0 J=\£-j\ 

, we can write the T-matrix element T^, e , n£ for an inelastic process 
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where P^ e (E) is the energy-dependent capture probability in the channel n for a given 
partial wave £ of a particular total angular momentum J, with similar definitions for 
P r J l , l f(E ) and P’l„ (E). Labeling the initial arrangement H 2 +D by a, the rate constant 
K n'^ n {E) = v ie i.a n ^ n (E) with u re i. = sj2E kin /= hk n /fj, a is simply 
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where f! — \ J — j '\,..., J+j' and t" = | J — j"\,..., J + j", respectively. As in the case of 
Eqs.Q and (jTlh we obtain the total quenching and reaction rate constant by summing 
over the appropriate channels and arrangements 
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(13) 

(14) 
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with P t J JE) = Pq(E) + JSE) + Er P„V(B). 

Due to the small mass of our system, we use the PST capture model to determine 
P£i i/{E) in the product channels only, where the kinetic energy is relatively large and 
the effect of near-threshold states below dynamical barriers may be neglected [52j, i.e., 


PAE) = { l] 


E — E n i > Vby 
otherwise 
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h 2 P(P + l )] 3 
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54^, Ce,, 

The van der Waals coefficient C 6n / and reduced mass /v appearing in the expression 
of the centrifugal barrier height V^y depends in the particular channel n' = ( a',v',j'). 
For PA(E) in the entrance channel, we use the low-energy expression 

PnAE) = A J n /E)kl M , 

where the coefficient _4,( f iE) is related to Jost. functions 


( 17 ) 
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Pj(2) 


At low-energy, J2e" Pn t"{P) is small and we can omit it in evaluating P Xot (E) ~ 
Pft(E). Furthermore, Pq, Pj{, and P t J ot are essentially ^-independent at low- 
energy; figures fl5l and fl6l show their variation with J for the initial state n = (v,j = 0) 
for para-H 2 and n = (v,j = 1) for ortho-H 2 , respectively. Although Pq, and P Xot 
all increase with J , the ratios Pn/Pt ot and Pr/ P^ ot are nearly constant (roughly 73% 
reaction and 27% quenching), allowing us to further simplify eq (jT3|) and (fT4l) as 


K Q(R)(£) = 



(18) 


where Pq^r)/Pua — PQ{R)/Ptot does not depend on J. In particular, their values 
for J = 0 is simply given by the corresponding number of open channels Nq for 
quenching and Nr for reaction, so that Pq(r)/P tot — Pm/PiLf = A'em/W... with 
iVtot = Nq + Nr. Although the rate constants for quenching and reaction processes 
depends on the details of the entrance channel via the sums over i and J, their ratio 
will not; the double-sum cancels out, and we find 


Kg<*>(E) P Q(R) Nq(r) 
K 9(E) Pq N q ' K «*(S) f,o, N m 


(19) 


Figured?] compares the results obtained using this statistical model and the full quantum 
treatment, We find that the agreement is roughly within 10-20% for v > 2, but quickly 
becomes worse when the reaction barrier plays an important role for very low v. Overall, 
this statistical model based on PST approximately gives the right branching ratio 
between quenching and reaction processes. 

As discussed in previous sections, the magnitude of the rate constants at low 
scattering energies are mainly given by the s-wave contribution: figures [5] and [6] show 
that this is the case for most of initial v up to 10 mK and higher. For para-H 2 in j = 0, 
the s-wave i = 0 implies J = 0 only, so that (2 J + 1)/(2j + 1) = 1, while for ortho-H 2 
in j = 1, £ = 0 implies J = 1 only, so that (2 J + 1)/(2j + 1) = 1 as well. In both cases, 
both sums over i and J in Eq. (fT8l) reduce to a single term; 

f A;{?£ 0 (A), for para-H 2 in j = 0 
x ’ . (20) 

[ A^jLo(E), for ortho-H 2 in j = 1 

For i = 0 with J = j, the coefficient N J n ^L {] (E — y 0) is simply related to the imaginary 
component j3 v j of the complex scattering length a v = a v ~ — ifi 
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channel n = ( a,v,j ) of the initial arrangement a=H 2 +D and total angular momentum 
J = j [55], 

A J n %„(E 40)= 4ft,, , (21) 

such that the rate constant takes the simple form 

J = j = 0, for para-H 2 

( 22 ) 

J — j — 1, for ortho-H 2 
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J J J J J 

Figure 15. Top panels: Pq(r) and Ptot f° r para-H 2 as a function of J for different 
entrance channels n = (j = 0,t>) with v = (v = 0 is not shown since 

only reaction is allowed, so that P^/P^ ot = 1). Lower panels: corresponding ratios 
PJ/PL- Also shown is the number of molecular states opened to quenching and 
reaction, Nq, and Nr respectively, and their ratios. These are J-independent and 
equal to the Pq^ values. 


Table [T] contains the values of a v j and f3 v j, which were extracted from the full quantum 
calculations [56] • They can be compared to an approximate value based on quantum 
reflection [57, 58] [59] , which is often used for barrierless reactions, 

1/4 
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(23) 


Using the value Cg -0 ’- 7-0 = 7.053 a.u. [44], together with fi a ~ 1836 a.u., we obtain 
/3 ~ 6.06 a.u., which is many order of magnitude larger than the values of /3 from the 
full quantum treatment. This illustrates well the fact that, although the statistical model 
gives approximately the right branching ratio of quenching to reaction processes, the 
approximation in Eq. (1251) cannot be used to obtain the magnitude of those processes. 
Although such an approximation can be applied to barrierless chemical systems, the 
reaction barrier reduces the value of j3 considerably, as shown in Table [D 
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Figure 17. Comparison of the branching ratios of quenching and reaction processes 
between the full quantum treatment Kn^ /(dash lines) and the statistical model 
NQ(R)/N tot (solid lines) for para-H 2 and ortho-H 2 . Results for v = 0 are not shown, 
since only reactive scattering is allowed. Ratios agree within 10-20%, except for v = 1 
for which the effect of the reaction barrier is more important. 


Figure 16. Same as Fig. [15] but for ortho-H 2 initially in j = 1. 
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Table 1. Real and imaginary contributions to the scattering length a v j = a v j — i(3 v ,j 
for para-H 2 (j = 0) and ortho-H 2 (j = 1) extracted from the quantum results. 


V 

D + para-H 2 (u, j = 0) 

D + ortho-H 2 (u, j = 1) 

a (a.u.) 

/ 3 (a.u.) 

a (a.u.) 

/ 3 (a.u.) 

0 

14.8 

7.9 x 10“ 14 

14.7 

2.8 x 10 -13 

1 

14.0 

2.4 x 10~ 7 

13.9 

3.1 x 10“ 7 

2 

13.4 

7.2 x 10 -5 

13.2 

1.1 x 10 -4 

3 

12.9 

1.2 x 10 -3 

12.7 

2.5 x 10 -3 

4 

12.2 

2.4 x 10 -2 

11.9 

7.3 x 10 -2 

5 

7.0 

3.0 x 10 _1 

6.6 

6.6 x 10" 1 


4. Conclusion 

The nuclear spin symmetry restricts the possible rotational states of H 2 , leading to para- 
142 with even j and ortho-H 2 with odd j. We investigated the effect of these restrictions 
on the scattering of H 2 with D, which can lead to elastic collisions, to quenching H 2 , 
or reaction to produce HD. We computed these various processes for the six lowest 
vibrational levels v of H 2 . We found structures in the energy range corresponding to 
100 mK and 10 K, due to low partial waves. Except for cases involving a sharp p -wave 
resonance, the results for para- and ortho-H 2 are very similar, with structures due to d 
or /-waves located at roughly the same energies in both cases, and with rate constants 
of comparable magnitude. 

There are marked differences for the cases in which a p -wave resonance is present, 
namely the three lowest v levels. For example, in the reactive case, while the sharp peak 
is present for both para- and ortho-H 2 in v — 0, it is absent for ortho-H 2 in v — 1 and 
para-H 2 in v = 2. The same is true for the quenching rate constant (though there is no 
quenching for v = 0, since it is the ground state). We also find the same overall behavior 
for the elastic cross sections. We also note that the exact position of the resonance is 
slightly affected by the initial value of j. 

We also investigated the state-to-state product formation, and found that the 
branching ratios are not sensitive to the scattering energy in the range considered 
here for para-H 2 . The same is true for ortho-H 2 except for the lower initial levels 
(v = 1, 2), where the p -wave resonance affects the branching ratios. By examining the 
rotationally-resolved rate constants of a specific initial case, namely v = 2, we can trace 
this dependence to sharper increases of the inelastic (either quenching or reaction) rate 
constant for certain rotational products near the resonance. These variations are due 
to the different coupling strengths of each total angular momentum J = 0, 1, and 2 
contributing to p -wave (£ = 1) initial partial wave for ortho-H 2 , as opposed to the case 
of para-H 2 , when only J = 1 contributes. 

The range of energies where the structure is present should be reachable 
experimentally, and thus could probe the effect of the nuclear symmetry on reaction 
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(or quenching). In addition, since the occurrence of p-wave resonances is very sensitive 
to the exact shape of the potential energy surface, detecting any structure in this 
benchmark system will provide a unique tool to compare with ab initio quantum 
chemistry calculations. Finally, based on our results, we suspect that the effect of the 
nuclear spin symmetry on cold and ultracold chemistry will be significant in systems 
displaying resonant features, and less so otherwise. However, if the nuclear spin is 
involved in the dynamics of the system, e.g., through an interaction term in the 
Hamiltonian, as in the studies of Refs. ggSSISZISSIHSlED], the nuclear spin symmetry 
could play an even more important role. Finally, we compared the results of our full 
quantum treatment to a statistical model, and found that although the model reproduces 
the proportion of quenching and reaction processes within 20% in most cases, it fails to 
give the magnitude of the rate constants. 
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